% From: http://www.ncbi.nlm.nih.gov/pubmed/15142845
% Title:Computer model of action potential of mouse ventricular myocytes
% Authors:Bondarenko, Szigeti, Bett, Kim, Rasmusson, 2004
function [VOI, STATES, ALGEBRAIC, CONSTANTS] = mainFunction()
% This is the "main function".  In Matlab, things work best if you rename this function to match the filename.
[VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel();
end

function [algebraicVariableCount] = getAlgebraicVariableCount()
% Used later when setting a global variable with the number of algebraic variables.
% Note: This is not the "main method".
algebraicVariableCount =71;
end
% There are a total of 41 entries in each of the rate and state variable arrays.
% There are a total of 73 entries in the constant variable array.


function [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel()
% Create ALGEBRAIC of correct size
global algebraicVariableCount;  algebraicVariableCount = getAlgebraicVariableCount();
% Initialise constants and state variables
[INIT_STATES, CONSTANTS] = initConsts;

% Set timespan to solve over
tspan = [0, 10];

% Set numerical accuracy options for ODE solver
options = odeset('RelTol', 1e-06, 'AbsTol', 1e-06, 'MaxStep', 1);

% Solve model with ODE solver
[VOI, STATES] = ode15s(@(VOI, STATES)computeRates(VOI, STATES, CONSTANTS), tspan, INIT_STATES, options);

% Compute algebraic variables
[RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS);
ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI);

% Plot state variables against variable of integration
[LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends();
figure();
plot(VOI, STATES);
xlabel(LEGEND_VOI);
l = legend(LEGEND_STATES);
set(l,'Interpreter','none');
end

function [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends()
LEGEND_STATES = ''; LEGEND_ALGEBRAIC = ''; LEGEND_VOI = ''; LEGEND_CONSTANTS = '';
LEGEND_t = strpad('time in component environment (millisecond)');
LEGEND_V_m = strpad('V in component membrane (millivolt)');
LEGEND_C_m = strpad('Cm in component membrane (microF_per_cm2)');
LEGEND_vol_myo = strpad('Vmyo in component membrane (microlitre)');
LEGEND_vol_JSR = strpad('VJSR in component membrane (microlitre)');
LEGEND_vol_NSR = strpad('VNSR in component membrane (microlitre)');
LEGEND_vol_ss = strpad('Vss in component membrane (microlitre)');
LEGEND_A_cap = strpad('Acap in component membrane (cm2)');
LEGEND_CONSTANTS(:,7) = strpad('Ko in component membrane (micromolar)');
LEGEND_Na_o = strpad('Nao in component membrane (micromolar)');
LEGEND_Ca_o = strpad('Cao in component membrane (micromolar)');
LEGEND_R = strpad('R in component membrane (joule_per_mole_kelvin)');
LEGEND_T = strpad('T in component membrane (kelvin)');
LEGEND_F = strpad('F in component membrane (coulomb_per_millimole)');
LEGEND_i_stim = strpad('i_stim in component membrane (picoA_per_picoF)');
LEGEND_i_CaL = strpad('i_CaL in component L_type_calcium_current (picoA_per_picoF)');
LEGEND_i_pCa = strpad('i_pCa in component calcium_pump_current (picoA_per_picoF)');
LEGEND_i_NaCa = strpad('i_NaCa in component sodium_calcium_exchange_current (picoA_per_picoF)');
LEGEND_i_CaB = strpad('i_Cab in component calcium_background_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,58) = strpad('i_Na in component fast_sodium_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,59) = strpad('i_Nab in component sodium_background_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,69) = strpad('i_NaK in component sodium_potassium_pump_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,61) = strpad('i_Kto_f in component fast_transient_outward_potassium_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,62) = strpad('i_Kto_s in component slow_transient_outward_potassium_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,63) = strpad('i_K1 in component time_independent_potassium_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,64) = strpad('i_Ks in component slow_delayed_rectifier_potassium_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,65) = strpad('i_Kur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,66) = strpad('i_Kss in component non_inactivating_steady_state_potassium_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,71) = strpad('i_ClCa in component calcium_activated_chloride_current (picoA_per_picoF)');
LEGEND_ALGEBRAIC(:,67) = strpad('i_Kr in component rapid_delayed_rectifier_potassium_current (picoA_per_picoF)');
LEGEND_CONSTANTS(:,13) = strpad('stim_start in component membrane (millisecond)');
LEGEND_stim_end = strpad('stim_end in component membrane (millisecond)');
LEGEND_stim_period = strpad('stim_period in component membrane (millisecond)');
LEGEND_stim_duration = strpad('stim_duration in component membrane (millisecond)');
LEGEND_stim_amp = strpad('stim_amplitude in component membrane (picoA_per_picoF)');
LEGEND_Ca_i = strpad('Cai in component calcium_concentration (micromolar)');
LEGEND_Ca_ss = strpad('Cass in component calcium_concentration (micromolar)');
LEGEND_Ca_JSR = strpad('CaJSR in component calcium_concentration (micromolar)');
LEGEND_Ca_NSR = strpad('CaNSR in component calcium_concentration (micromolar)');
LEGEND_B_i = strpad('Bi in component calcium_concentration (dimensionless)');
LEGEND_B_ss = strpad('Bss in component calcium_concentration (dimensionless)');
LEGEND_B_JSR = strpad('BJSR in component calcium_concentration (dimensionless)');
LEGEND_CMDN_tot = strpad('CMDN_tot in component calcium_concentration (micromolar)');
LEGEND_CSQN_tot = strpad('CSQN_tot in component calcium_concentration (micromolar)');
LEGEND_Km_CMDN = strpad('Km_CMDN in component calcium_concentration (micromolar)');
LEGEND_CONSTANTS(:,21) = strpad('Km_CSQN in component calcium_concentration (micromolar)');
LEGEND_J_leak = strpad('J_leak in component calcium_fluxes (micromolar_per_millisecond)');
LEGEND_J_rel = strpad('J_rel in component calcium_fluxes (micromolar_per_millisecond)');
LEGEND_J_up = strpad('J_up in component calcium_fluxes (micromolar_per_millisecond)');
LEGEND_J_tr = strpad('J_tr in component calcium_fluxes (micromolar_per_millisecond)');
LEGEND_J_trpn = strpad('J_trpn in component calcium_fluxes (micromolar_per_millisecond)');
LEGEND_J_xfer = strpad('J_xfer in component calcium_fluxes (micromolar_per_millisecond)');
LEGEND_CONSTANTS(:,22) = strpad('k_plus_htrpn in component calcium_fluxes (per_micromolar_millisecond)');
LEGEND_CONSTANTS(:,23) = strpad('k_minus_htrpn in component calcium_fluxes (per_millisecond)');
LEGEND_CONSTANTS(:,24) = strpad('k_plus_ltrpn in component calcium_fluxes (per_micromolar_millisecond)');
LEGEND_CONSTANTS(:,25) = strpad('k_minus_ltrpn in component calcium_fluxes (per_millisecond)');
LEGEND_P_RyR = strpad('P_RyR in component calcium_fluxes (dimensionless)');
LEGEND_v1 = strpad('v1 in component calcium_fluxes (per_millisecond)');
LEGEND_tau_tr = strpad('tau_tr in component calcium_fluxes (millisecond)');
LEGEND_v2 = strpad('v2 in component calcium_fluxes (per_millisecond)');
LEGEND_tau_xfer = strpad('tau_xfer in component calcium_fluxes (millisecond)');
LEGEND_v3 = strpad('v3 in component calcium_fluxes (micromolar_per_millisecond)');
LEGEND_K_m_up = strpad('Km_up in component calcium_fluxes (micromolar)');
LEGEND_CONSTANTS(:,32) = strpad('LTRPN_tot in component calcium_buffering (micromolar)');
LEGEND_CONSTANTS(:,33) = strpad('HTRPN_tot in component calcium_buffering (micromolar)');
LEGEND_STATES(:,7) = strpad('LTRPN_Ca in component calcium_buffering (micromolar)');
LEGEND_STATES(:,8) = strpad('HTRPN_Ca in component calcium_buffering (micromolar)');
LEGEND_i_CaL_max= strpad('i_CaL_max in component L_type_calcium_current (picoA_per_picoF)');
LEGEND_P_O1 = strpad('P_O1 in component ryanodine_receptors (dimensionless)');
LEGEND_P_O2 = strpad('P_O2 in component ryanodine_receptors (dimensionless)');
LEGEND_P_C1 = strpad('P_C1 in component ryanodine_receptors (dimensionless)');
LEGEND_P_C2 = strpad('P_C2 in component ryanodine_receptors (dimensionless)');
LEGEND_k_plus_a = strpad('k_plus_a in component ryanodine_receptors (micromolar4_per_millisecond)');
LEGEND_k_minus_a = strpad('k_minus_a in component ryanodine_receptors (per_millisecond)');
LEGEND_k_plus_b = strpad('k_plus_b in component ryanodine_receptors (micromolar3_per_millisecond)');
LEGEND_k_minus_b = strpad('k_minus_b in component ryanodine_receptors (per_millisecond)');
LEGEND_k_plus_c = strpad('k_plus_c in component ryanodine_receptors (per_millisecond)');
LEGEND_k_minus_c = strpad('k_minus_c in component ryanodine_receptors (per_millisecond)');
LEGEND_m = strpad('m in component ryanodine_receptors (dimensionless)');
LEGEND_n = strpad('n in component ryanodine_receptors (dimensionless)');
LEGEND_E_CaL = strpad('E_CaL in component L_type_calcium_current (millivolt)');
LEGEND_G_CaL = strpad('g_CaL in component L_type_calcium_current (milliS_per_microF)');
LEGEND_O = strpad('O in component L_type_calcium_current (dimensionless)');
LEGEND_C1 = strpad('C1 in component L_type_calcium_current (dimensionless)');
LEGEND_C2 = strpad('C2 in component L_type_calcium_current (dimensionless)');
LEGEND_C3 = strpad('C3 in component L_type_calcium_current (dimensionless)');
LEGEND_C4 = strpad('C4 in component L_type_calcium_current (dimensionless)');
LEGEND_I1 = strpad('I1 in component L_type_calcium_current (dimensionless)');
LEGEND_I2 = strpad('I2 in component L_type_calcium_current (dimensionless)');
LEGEND_I3 = strpad('I3 in component L_type_calcium_current (dimensionless)');
LEGEND_alpha = strpad('alpha in component L_type_calcium_current (per_millisecond)');
LEGEND_beta = strpad('beta in component L_type_calcium_current (per_millisecond)');
LEGEND_gamma = strpad('gamma in component L_type_calcium_current (per_millisecond)');
LEGEND_K_pcf = strpad('Kpcf in component L_type_calcium_current (per_millisecond)');
LEGEND_K_pcb = strpad('Kpcb in component L_type_calcium_current (per_millisecond)');
LEGEND_k_pc_max = strpad('Kpc_max in component L_type_calcium_current (per_millisecond)');
LEGEND_k_pc_half = strpad('Kpc_half in component L_type_calcium_current (micromolar)');
LEGEND_i_pCa_max = strpad('i_pCa_max in component calcium_pump_current (picoA_per_picoF)');
LEGEND_Km_pCa= strpad('Km_pCa in component calcium_pump_current (micromolar)');
LEGEND_k_NaCa = strpad('k_NaCa in component sodium_calcium_exchange_current (picoA_per_picoF)');
LEGEND_K_mNa = strpad('K_mNa in component sodium_calcium_exchange_current (micromolar)');
LEGEND_K_mCa = strpad('K_mCa in component sodium_calcium_exchange_current (micromolar)');
LEGEND_k_sat = strpad('k_sat in component sodium_calcium_exchange_current (dimensionless)');
LEGEND_eta = strpad('eta in component sodium_calcium_exchange_current (dimensionless)');
LEGEND_Na_i = strpad('Nai in component sodium_concentration (micromolar)');
LEGEND_G_CaB = strpad('g_Cab in component calcium_background_current (milliS_per_microF)');
LEGEND_E_CaN = strpad('E_CaN in component calcium_background_current (millivolt)');
LEGEND_ALGEBRAIC(:,57) = strpad('E_Na in component fast_sodium_current (millivolt)');
LEGEND_CONSTANTS(:,56) = strpad('g_Na in component fast_sodium_current (milliS_per_microF)');
LEGEND_STATES(:,20) = strpad('O_Na in component fast_sodium_current (dimensionless)');
LEGEND_STATES(:,21) = strpad('C_Na1 in component fast_sodium_current (dimensionless)');
LEGEND_STATES(:,22) = strpad('C_Na2 in component fast_sodium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,4) = strpad('C_Na3 in component fast_sodium_current (dimensionless)');
LEGEND_STATES(:,23) = strpad('I1_Na in component fast_sodium_current (dimensionless)');
LEGEND_STATES(:,24) = strpad('I2_Na in component fast_sodium_current (dimensionless)');
LEGEND_STATES(:,25) = strpad('IF_Na in component fast_sodium_current (dimensionless)');
LEGEND_STATES(:,26) = strpad('IC_Na2 in component fast_sodium_current (dimensionless)');
LEGEND_STATES(:,27) = strpad('IC_Na3 in component fast_sodium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,14) = strpad('alpha_Na11 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,36) = strpad('beta_Na11 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,27) = strpad('alpha_Na12 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,38) = strpad('beta_Na12 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,32) = strpad('alpha_Na13 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,40) = strpad('beta_Na13 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,42) = strpad('alpha_Na3 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,44) = strpad('beta_Na3 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,46) = strpad('alpha_Na2 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,48) = strpad('beta_Na2 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,50) = strpad('alpha_Na4 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,52) = strpad('beta_Na4 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,54) = strpad('alpha_Na5 in component fast_sodium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,56) = strpad('beta_Na5 in component fast_sodium_current (per_millisecond)');
LEGEND_STATES(:,28) = strpad('Ki in component potassium_concentration (micromolar)');
LEGEND_CONSTANTS(:,57) = strpad('g_Nab in component sodium_background_current (milliS_per_microF)');
LEGEND_ALGEBRAIC(:,60) = strpad('E_K in component fast_transient_outward_potassium_current (millivolt)');
LEGEND_CONSTANTS(:,58) = strpad('g_Kto_f in component fast_transient_outward_potassium_current (milliS_per_microF)');
LEGEND_STATES(:,29) = strpad('ato_f in component fast_transient_outward_potassium_current (dimensionless)');
LEGEND_STATES(:,30) = strpad('ito_f in component fast_transient_outward_potassium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,5) = strpad('alpha_a in component fast_transient_outward_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,15) = strpad('beta_a in component fast_transient_outward_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,6) = strpad('alpha_i in component fast_transient_outward_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,16) = strpad('beta_i in component fast_transient_outward_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,7) = strpad('ass in component slow_transient_outward_potassium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,8) = strpad('iss in component slow_transient_outward_potassium_current (dimensionless)');
LEGEND_CONSTANTS(:,59) = strpad('g_Kto_s in component slow_transient_outward_potassium_current (milliS_per_microF)');
LEGEND_STATES(:,31) = strpad('ato_s in component slow_transient_outward_potassium_current (dimensionless)');
LEGEND_STATES(:,32) = strpad('ito_s in component slow_transient_outward_potassium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,17) = strpad('tau_ta_s in component slow_transient_outward_potassium_current (millisecond)');
LEGEND_ALGEBRAIC(:,18) = strpad('tau_ti_s in component slow_transient_outward_potassium_current (millisecond)');
LEGEND_CONSTANTS(:,60) = strpad('g_Ks in component slow_delayed_rectifier_potassium_current (milliS_per_microF)');
LEGEND_STATES(:,33) = strpad('nKs in component slow_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,9) = strpad('alpha_n in component slow_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,19) = strpad('beta_n in component slow_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_CONSTANTS(:,61) = strpad('g_Kur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (milliS_per_microF)');
LEGEND_STATES(:,34) = strpad('aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_STATES(:,35) = strpad('iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,20) = strpad('tau_aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (millisecond)');
LEGEND_ALGEBRAIC(:,21) = strpad('tau_iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (millisecond)');
LEGEND_CONSTANTS(:,62) = strpad('g_Kss in component non_inactivating_steady_state_potassium_current (milliS_per_microF)');
LEGEND_STATES(:,36) = strpad('aKss in component non_inactivating_steady_state_potassium_current (dimensionless)');
LEGEND_STATES(:,37) = strpad('iKss in component non_inactivating_steady_state_potassium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,22) = strpad('tau_Kss in component non_inactivating_steady_state_potassium_current (millisecond)');
LEGEND_CONSTANTS(:,63) = strpad('g_Kr in component rapid_delayed_rectifier_potassium_current (milliS_per_microF)');
LEGEND_STATES(:,38) = strpad('O_K in component rapid_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_STATES(:,39) = strpad('C_K1 in component rapid_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_STATES(:,40) = strpad('C_K2 in component rapid_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,10) = strpad('C_K0 in component rapid_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_STATES(:,41) = strpad('I_K in component rapid_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_ALGEBRAIC(:,23) = strpad('alpha_a0 in component rapid_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,28) = strpad('beta_a0 in component rapid_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_CONSTANTS(:,64) = strpad('kb in component rapid_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_CONSTANTS(:,65) = strpad('kf in component rapid_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,11) = strpad('alpha_a1 in component rapid_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,24) = strpad('beta_a1 in component rapid_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,29) = strpad('alpha_i in component rapid_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_ALGEBRAIC(:,33) = strpad('beta_i in component rapid_delayed_rectifier_potassium_current (per_millisecond)');
LEGEND_CONSTANTS(:,66) = strpad('i_NaK_max in component sodium_potassium_pump_current (picoA_per_picoF)');
LEGEND_CONSTANTS(:,67) = strpad('Km_Nai in component sodium_potassium_pump_current (micromolar)');
LEGEND_CONSTANTS(:,68) = strpad('Km_Ko in component sodium_potassium_pump_current (micromolar)');
LEGEND_ALGEBRAIC(:,68) = strpad('f_NaK in component sodium_potassium_pump_current (dimensionless)');
LEGEND_CONSTANTS(:,72) = strpad('sigma in component sodium_potassium_pump_current (dimensionless)');
LEGEND_CONSTANTS(:,69) = strpad('g_ClCa in component calcium_activated_chloride_current (milliS_per_microF)');
LEGEND_ALGEBRAIC(:,70) = strpad('O_ClCa in component calcium_activated_chloride_current (dimensionless)');
LEGEND_CONSTANTS(:,70) = strpad('E_Cl in component calcium_activated_chloride_current (millivolt)');
LEGEND_CONSTANTS(:,71) = strpad('Km_Cl in component calcium_activated_chloride_current (micromolar)');
LEGEND_V_m_prime = strpad('d/dt V in component membrane (millivolt)');
LEGEND_Ca_i_prime = strpad('d/dt Cai in component calcium_concentration (micromolar)');
LEGEND_Ca_ss_prime = strpad('d/dt Cass in component calcium_concentration (micromolar)');
LEGEND_Ca_JSR_prime = strpad('d/dt CaJSR in component calcium_concentration (micromolar)');
LEGEND_Ca_NSR_prime = strpad('d/dt CaNSR in component calcium_concentration (micromolar)');
LEGEND_P_RyR_prime = strpad('d/dt P_RyR in component calcium_fluxes (dimensionless)');
LEGEND_RATES(:,7) = strpad('d/dt LTRPN_Ca in component calcium_buffering (micromolar)');
LEGEND_RATES(:,8) = strpad('d/dt HTRPN_Ca in component calcium_buffering (micromolar)');
LEGEND_P_O1_prime = strpad('d/dt P_O1 in component ryanodine_receptors (dimensionless)');
LEGEND_P_O2_prime = strpad('d/dt P_O2 in component ryanodine_receptors (dimensionless)');
LEGEND_P_C2_prime = strpad('d/dt P_C2 in component ryanodine_receptors (dimensionless)');
LEGEND_O_prime = strpad('d/dt O in component L_type_calcium_current (dimensionless)');
LEGEND_C2_prime = strpad('d/dt C2 in component L_type_calcium_current (dimensionless)');
LEGEND_C3_prime = strpad('d/dt C3 in component L_type_calcium_current (dimensionless)');
LEGEND_C4_prime = strpad('d/dt C4 in component L_type_calcium_current (dimensionless)');
LEGEND_I1_prime = strpad('d/dt I1 in component L_type_calcium_current (dimensionless)');
LEGEND_I2_prime = strpad('d/dt I2 in component L_type_calcium_current (dimensionless)');
LEGEND_I3_prime = strpad('d/dt I3 in component L_type_calcium_current (dimensionless)');
LEGEND_RATES(:,19) = strpad('d/dt Nai in component sodium_concentration (micromolar)');
LEGEND_RATES(:,22) = strpad('d/dt C_Na2 in component fast_sodium_current (dimensionless)');
LEGEND_RATES(:,21) = strpad('d/dt C_Na1 in component fast_sodium_current (dimensionless)');
LEGEND_RATES(:,20) = strpad('d/dt O_Na in component fast_sodium_current (dimensionless)');
LEGEND_RATES(:,25) = strpad('d/dt IF_Na in component fast_sodium_current (dimensionless)');
LEGEND_RATES(:,23) = strpad('d/dt I1_Na in component fast_sodium_current (dimensionless)');
LEGEND_RATES(:,24) = strpad('d/dt I2_Na in component fast_sodium_current (dimensionless)');
LEGEND_RATES(:,26) = strpad('d/dt IC_Na2 in component fast_sodium_current (dimensionless)');
LEGEND_RATES(:,27) = strpad('d/dt IC_Na3 in component fast_sodium_current (dimensionless)');
LEGEND_RATES(:,28) = strpad('d/dt Ki in component potassium_concentration (micromolar)');
LEGEND_RATES(:,29) = strpad('d/dt ato_f in component fast_transient_outward_potassium_current (dimensionless)');
LEGEND_RATES(:,30) = strpad('d/dt ito_f in component fast_transient_outward_potassium_current (dimensionless)');
LEGEND_RATES(:,31) = strpad('d/dt ato_s in component slow_transient_outward_potassium_current (dimensionless)');
LEGEND_RATES(:,32) = strpad('d/dt ito_s in component slow_transient_outward_potassium_current (dimensionless)');
LEGEND_RATES(:,33) = strpad('d/dt nKs in component slow_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_RATES(:,34) = strpad('d/dt aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_RATES(:,35) = strpad('d/dt iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_RATES(:,36) = strpad('d/dt aKss in component non_inactivating_steady_state_potassium_current (dimensionless)');
LEGEND_RATES(:,37) = strpad('d/dt iKss in component non_inactivating_steady_state_potassium_current (dimensionless)');
LEGEND_RATES(:,40) = strpad('d/dt C_K2 in component rapid_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_RATES(:,39) = strpad('d/dt C_K1 in component rapid_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_RATES(:,38) = strpad('d/dt O_K in component rapid_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_RATES(:,41) = strpad('d/dt I_K in component rapid_delayed_rectifier_potassium_current (dimensionless)');
LEGEND_STATES  = LEGEND_STATES';
LEGEND_ALGEBRAIC = LEGEND_ALGEBRAIC';
LEGEND_RATES = LEGEND_RATES';
LEGEND_CONSTANTS = LEGEND_CONSTANTS';
end

function [STATES, CONSTANTS] = initConsts()
t = 0; CONSTANTS = []; STATES = []; ALGEBRAIC = [];
V_m = -82.4202;
C_m = 1;
vol_myo = 25.84e-6;
vol_JSR = 0.12e-6;
vol_NSR = 2.098e-6;
vol_ss = 1.485e-9;
A_cap = 1.534e-4;
CONSTANTS(:,7) = 5400;
Na_o = 140000;
Ca_o = 1800;
R = 8.314;
T = 298;
F = 96.5;
stim_start = 20;
stim_end = 100000;
stim_period = 71.43;
stim_duration = 0.5;
stim_amp = -80;
Ca_i = 0.115001;
Ca_ss = 0.115001;
Ca_JSR = 1299.5;
Ca_NSR = 1299.5;
CMDN_tot = 50;
CSQN_tot = 15000;
Km_CMDN = 0.238;
Km_CSQN = 800;
CONSTANTS(:,22) = 0.00237;
CONSTANTS(:,23) = 3.2e-5;
CONSTANTS(:,24) = 0.0327;
CONSTANTS(:,25) = 0.0196;
P_RyR = 0;
v1 = 4.5;
tau_tr = 20;
v2 = 1.74e-5;
tau_xfer = 8;
v3 = 0.45;
K_m_up = 0.5;
CONSTANTS(:,32) = 70;
CONSTANTS(:,33) = 140;
STATES(:,7) = 11.2684;
STATES(:,8) = 125.29;
i_CaL_max= 7;
P_O1 = 0.149102e-4;
P_O2 = 0.951726e-10;
P_C2 = 0.16774e-3;
k_plus_a = 0.006075;
k_minus_a = 0.07125;
k_plus_b = 0.00405;
k_minus_b = 0.965;
k_plus_c = 0.009;
k_minus_c = 0.0008;
m = 3;
n = 4;
E_CaL = 63;
G_CaL = 0.1729;
O = 0.930308e-18;
C2 = 0.124216e-3;
C3 = 0.578679e-8;
C4 = 0.119816e-12;
I1 = 0.497923e-18;
I2 = 0.345847e-13;
I3 = 0.185106e-13;
K_pcb = 0.0005;
k_pc_max = 0.23324;
k_pc_half = 20;
i_pCa_max = 1;
Km_pCa= 0.5;
k_NaCa = 292.8;
K_mNa = 87500;
K_mCa = 1380;
k_sat = 0.1;
eta = 0.35;
Na_i = 14237.1;
G_CaB = 0.000367;
CONSTANTS(:,56) = 13;
STATES(:,20) = 0.713483e-6;
STATES(:,21) = 0.279132e-3;
STATES(:,22) = 0.020752;
STATES(:,23) = 0.673345e-6;
STATES(:,24) = 0.155787e-8;
STATES(:,25) = 0.153176e-3;
STATES(:,26) = 0.0113879;
STATES(:,27) = 0.34278;
STATES(:,28) = 143720;
CONSTANTS(:,57) = 0.0026;
CONSTANTS(:,58) = 0.4067;
STATES(:,29) = 0.265563e-2;
STATES(:,30) = 0.999977;
CONSTANTS(:,59) = 0;
STATES(:,31) = 0.417069e-3;
STATES(:,32) = 0.998543;
CONSTANTS(:,60) = 0.00575;
STATES(:,33) = 0.262753e-3;
CONSTANTS(:,61) = 0.16;
STATES(:,34) = 0.417069e-3;
STATES(:,35) = 0.998543;
CONSTANTS(:,62) = 0.05;
STATES(:,36) = 0.417069e-3;
STATES(:,37) = 1;
CONSTANTS(:,63) = 0.078;
STATES(:,38) = 0.175298e-3;
STATES(:,39) = 0.992513e-3;
STATES(:,40) = 0.641229e-3;
STATES(:,41) = 0.319129e-4;
CONSTANTS(:,64) = 0.036778;
CONSTANTS(:,65) = 0.023761;
CONSTANTS(:,66) = 0.88;
CONSTANTS(:,67) = 21000;
CONSTANTS(:,68) = 1500;
CONSTANTS(:,69) = 10;
CONSTANTS(:,70) = -40;
CONSTANTS(:,71) = 10;
CONSTANTS(:,72) =  (1.00000./7.00000).*((exp((Na_o./67300.0))) - 1.00000);
CONSTANTS(:,72) = 0.00000;
if (isempty(STATES)), warning('Initial values for states not set');, end
end

function [RATES, ALGEBRAIC] = computeRates(t, STATES, CONSTANTS)
global algebraicVariableCount;
statesSize = size(STATES);
statesColumnCount = statesSize(2);
if ( statesColumnCount == 1)
    STATES = STATES';
    ALGEBRAIC = zeros(1, algebraicVariableCount);
else
    statesRowCount = statesSize(1);
    ALGEBRAIC = zeros(statesRowCount, algebraicVariableCount);
    RATES = zeros(statesRowCount, statesColumnCount);
end
RATES(:,37) = CONSTANTS(:,72);
RATES(:,7) =  CONSTANTS(:,24).*Ca_i.*(CONSTANTS(:,32) - STATES(:,7)) -  CONSTANTS(:,25).*STATES(:,7);
RATES(:,8) =  CONSTANTS(:,22).*Ca_i.*(CONSTANTS(:,33) - STATES(:,8)) -  CONSTANTS(:,23).*STATES(:,8);
P_O2_prime =  k_plus_b.*(Ca_ss .^ m).*P_O1 -  k_minus_b.*P_O2;
P_C2_prime =  k_plus_c.*P_O1 -  k_minus_c.*P_C2;
P_C1 = 1.00000 - (P_C2+P_O1+P_O2);
P_O1_prime = ( k_plus_a.*(Ca_ss .^ n).*P_C1+ k_minus_b.*P_O2+ k_minus_c.*P_C2) - ( k_minus_a.*P_O1+ k_plus_b.*(Ca_ss .^ m).*P_O1+ k_plus_c.*P_O1);
ALGEBRAIC(:,5) =  0.180640.*(exp(( 0.0357700.*(V_m+30.0000))));
ALGEBRAIC(:,15) =  0.395600.*(exp((  - 0.0623700.*(V_m+30.0000))));
RATES(:,29) =  ALGEBRAIC(:,5).*(1.00000 - STATES(:,29)) -  ALGEBRAIC(:,15).*STATES(:,29);
ALGEBRAIC(:,6) = ( 0.000152000.*(exp(( - (V_m+13.5000)./7.00000))))./( 0.00670830.*(exp(( - (V_m+33.5000)./7.00000)))+1.00000);
ALGEBRAIC(:,16) = ( 0.000950000.*(exp(((V_m+33.5000)./7.00000))))./( 0.0513350.*(exp(((V_m+33.5000)./7.00000)))+1.00000);
RATES(:,30) =  ALGEBRAIC(:,6).*(1.00000 - STATES(:,30)) -  ALGEBRAIC(:,16).*STATES(:,30);
ALGEBRAIC(:,7) = 1.00000./(1.00000+(exp(( - (V_m+22.5000)./7.70000))));
ALGEBRAIC(:,17) =  0.493000.*(exp((  - 0.0629000.*V_m)))+2.05800;
RATES(:,31) = (ALGEBRAIC(:,7) - STATES(:,31))./ALGEBRAIC(:,17);
ALGEBRAIC(:,8) = 1.00000./(1.00000+(exp(((V_m+45.2000)./5.70000))));
ALGEBRAIC(:,18) = 270.000+1050.00./(1.00000+(exp(((V_m+45.2000)./5.70000))));
RATES(:,32) = (ALGEBRAIC(:,8) - STATES(:,32))./ALGEBRAIC(:,18);
ALGEBRAIC(:,9) = ( 4.81333e-06.*(V_m+26.5000))./(1.00000 - (exp((  - 0.128000.*(V_m+26.5000)))));
ALGEBRAIC(:,19) =  9.53333e-05.*(exp((  - 0.0380000.*(V_m+26.5000))));
RATES(:,33) =  ALGEBRAIC(:,9).*(1.00000 - STATES(:,33)) -  ALGEBRAIC(:,19).*STATES(:,33);
ALGEBRAIC(:,20) =  0.493000.*(exp((  - 0.0629000.*V_m)))+2.05800;
RATES(:,34) = (ALGEBRAIC(:,7) - STATES(:,34))./ALGEBRAIC(:,20);
ALGEBRAIC(:,21) = 1200.00 - 170.000./(1.00000+(exp(((V_m+45.2000)./5.70000))));
RATES(:,35) = (ALGEBRAIC(:,8) - STATES(:,35))./ALGEBRAIC(:,21);
ALGEBRAIC(:,22) =  39.3000.*(exp((  - 0.0862000.*V_m)))+13.1700;
RATES(:,36) = (ALGEBRAIC(:,7) - STATES(:,36))./ALGEBRAIC(:,22);
ALGEBRAIC(:,11) =  0.0137330.*(exp(( 0.0381980.*V_m)));
ALGEBRAIC(:,24) =  6.89000e-05.*(exp((  - 0.0417800.*V_m)));
RATES(:,40) = ( CONSTANTS(:,65).*STATES(:,39)+ ALGEBRAIC(:,24).*STATES(:,38)) - ( CONSTANTS(:,64).*STATES(:,40)+ ALGEBRAIC(:,11).*STATES(:,40));
C1 = 1.00000 - (O+C2+C3+C4+I1+I2+I3);
alpha = ( 0.400000.*(exp(((V_m+12.0000)./10.0000))).*((1.00000+ 0.700000.*(exp(( - ((V_m+40.0000) .^ 2.00000)./10.0000)))) -  0.750000.*(exp(( - ((V_m+20.0000) .^ 2.00000)./400.000)))))./(1.00000+ 0.120000.*(exp(((V_m+12.0000)./10.0000))));
beta =  0.0500000.*(exp(( - (V_m+12.0000)./13.0000)));
C2_prime = ( 4.00000.*alpha.*C1+ 2.00000.*beta.*C3) - ( beta.*C2+ 3.00000.*alpha.*C2);
C3_prime = ( 3.00000.*alpha.*C2+ 3.00000.*beta.*C4) - ( 2.00000.*beta.*C3+ 2.00000.*alpha.*C3);
ALGEBRAIC(:,10) = 1.00000 - (STATES(:,39)+STATES(:,40)+STATES(:,38)+STATES(:,41));
ALGEBRAIC(:,23) =  0.0223480.*(exp(( 0.0117600.*V_m)));
ALGEBRAIC(:,28) =  0.0470020.*(exp((  - 0.0631000.*V_m)));
RATES(:,39) = ( ALGEBRAIC(:,23).*ALGEBRAIC(:,10)+ CONSTANTS(:,64).*STATES(:,40)) - ( ALGEBRAIC(:,28).*STATES(:,39)+ CONSTANTS(:,65).*STATES(:,39));
ALGEBRAIC(:,29) =  0.0908210.*(exp(( 0.0233910.*(V_m+5.00000))));
ALGEBRAIC(:,33) =  0.00649700.*(exp((  - 0.0326800.*(V_m+5.00000))));
RATES(:,38) = ( ALGEBRAIC(:,11).*STATES(:,40)+ ALGEBRAIC(:,33).*STATES(:,41)) - ( ALGEBRAIC(:,24).*STATES(:,38)+ ALGEBRAIC(:,29).*STATES(:,38));
RATES(:,41) =  ALGEBRAIC(:,29).*STATES(:,38) -  ALGEBRAIC(:,33).*STATES(:,41);
gamma = ( k_pc_max.*Ca_ss)./(k_pc_half+Ca_ss);
K_pcf =  13.0000.*(1.00000 - (exp(( - ((V_m+14.5000) .^ 2.00000)./100.000))));
O_prime = ( alpha.*C4+ K_pcb.*I1+ 0.00100000.*( alpha.*I2 -  K_pcf.*O)) - ( 4.00000.*beta.*O+ gamma.*O);
C4_prime = ( 2.00000.*alpha.*C3+ 4.00000.*beta.*O+ 0.0100000.*( 4.00000.*K_pcb.*beta.*I1 -  alpha.*gamma.*C4)+ 0.00200000.*( 4.00000.*beta.*I2 -  K_pcf.*C4)+ 4.00000.*beta.*K_pcb.*I3) - ( 3.00000.*beta.*C4+ alpha.*C4+ 1.00000.*gamma.*K_pcf.*C4);
I1_prime = ( gamma.*O+ 0.00100000.*( alpha.*I3 -  K_pcf.*I1)+ 0.0100000.*( alpha.*gamma.*C4 -  4.00000.*beta.*K_pcf.*I1)) -  K_pcb.*I1;
I2_prime = ( 0.00100000.*( K_pcf.*O -  alpha.*I2)+ K_pcb.*I3+ 0.00200000.*( K_pcf.*C4 -  4.00000.*beta.*I2)) -  gamma.*I2;
I3_prime = ( 0.00100000.*( K_pcf.*I1 -  alpha.*I3)+ gamma.*I2+ 1.00000.*gamma.*K_pcf.*C4) - ( 4.00000.*beta.*K_pcb.*I3+ K_pcb.*I3);
B_JSR = (1.00000+( CSQN_tot.*Km_CSQN)./((Km_CSQN+Ca_JSR) .^ 2.00000)) .^ ( - 1.00000);
J_rel =  v1.*(P_O1+P_O2).*(Ca_JSR - Ca_ss).*P_RyR;
J_tr = (Ca_NSR - Ca_JSR)./tau_tr;
Ca_JSR_prime =  B_JSR.*(J_tr - J_rel);
J_leak =  v2.*(Ca_NSR - Ca_i);
J_up = ( v3.*(Ca_i .^ 2.00000))./((K_m_up .^ 2.00000)+(Ca_i .^ 2.00000));
Ca_NSR_prime = ( (J_up - J_leak).*vol_myo)./vol_NSR - ( J_tr.*vol_JSR)./vol_NSR;
ALGEBRAIC(:,4) = 1.00000 - (STATES(:,20)+STATES(:,21)+STATES(:,22)+STATES(:,25)+STATES(:,23)+STATES(:,24)+STATES(:,26)+STATES(:,27));
ALGEBRAIC(:,14) = 3.80200./( 0.102700.*(exp(( - (V_m+2.50000)./17.0000)))+ 0.200000.*(exp(( - (V_m+2.50000)./150.000))));
ALGEBRAIC(:,36) =  0.191700.*(exp(( - (V_m+2.50000)./20.3000)));
ALGEBRAIC(:,27) = 3.80200./( 0.102700.*(exp(( - (V_m+2.50000)./15.0000)))+ 0.230000.*(exp(( - (V_m+2.50000)./150.000))));
ALGEBRAIC(:,38) =  0.200000.*(exp(( - (V_m - 2.50000)./20.3000)));
ALGEBRAIC(:,42) =  7.00000e-07.*(exp(( - (V_m+7.00000)./7.70000)));
ALGEBRAIC(:,44) = 0.00840000+ 2.00000e-05.*(V_m+7.00000);
RATES(:,22) = ( ALGEBRAIC(:,14).*ALGEBRAIC(:,4)+ ALGEBRAIC(:,38).*STATES(:,21)+ ALGEBRAIC(:,42).*STATES(:,26)) - ( ALGEBRAIC(:,36).*STATES(:,22)+ ALGEBRAIC(:,27).*STATES(:,22)+ ALGEBRAIC(:,44).*STATES(:,22));
ALGEBRAIC(:,32) = 3.80200./( 0.102700.*(exp(( - (V_m+2.50000)./12.0000)))+ 0.250000.*(exp(( - (V_m+2.50000)./150.000))));
ALGEBRAIC(:,40) =  0.220000.*(exp(( - (V_m - 7.50000)./20.3000)));
RATES(:,21) = ( ALGEBRAIC(:,27).*STATES(:,22)+ ALGEBRAIC(:,40).*STATES(:,20)+ ALGEBRAIC(:,42).*STATES(:,25)) - ( ALGEBRAIC(:,38).*STATES(:,21)+ ALGEBRAIC(:,32).*STATES(:,21)+ ALGEBRAIC(:,44).*STATES(:,21));
RATES(:,26) = ( ALGEBRAIC(:,14).*STATES(:,27)+ ALGEBRAIC(:,38).*STATES(:,25)+ ALGEBRAIC(:,44).*STATES(:,22)) - ( ALGEBRAIC(:,36).*STATES(:,26)+ ALGEBRAIC(:,27).*STATES(:,26)+ ALGEBRAIC(:,42).*STATES(:,26));
RATES(:,27) = ( ALGEBRAIC(:,36).*STATES(:,26)+ ALGEBRAIC(:,44).*ALGEBRAIC(:,4)) - ( ALGEBRAIC(:,14).*STATES(:,27)+ ALGEBRAIC(:,42).*STATES(:,27));
i_CaL =  G_CaL.*O.*(V_m - E_CaL);
B_ss = (1.00000+( CMDN_tot.*Km_CMDN)./((Km_CMDN+Ca_ss) .^ 2.00000)) .^ ( - 1.00000);
J_xfer = (Ca_ss - Ca_i)./tau_xfer;
Ca_ss_prime =  B_ss.*(( J_rel.*vol_JSR)./vol_ss - (( J_xfer.*vol_myo)./vol_ss+( i_CaL.*A_cap.*C_m)./( 2.00000.*vol_ss.*F)));
P_RyR_prime =   - 0.0400000.*P_RyR -  (( 0.100000.*i_CaL)./i_CaL_max).*(exp(( - ((V_m - 5.00000) .^ 2.00000)./648.000)));
ALGEBRAIC(:,46) = 1.00000./( 0.188495.*(exp(( - (V_m+7.00000)./16.6000)))+0.393956);
ALGEBRAIC(:,48) = ( ALGEBRAIC(:,32).*ALGEBRAIC(:,46).*ALGEBRAIC(:,42))./( ALGEBRAIC(:,40).*ALGEBRAIC(:,44));
RATES(:,20) = ( ALGEBRAIC(:,32).*STATES(:,21)+ ALGEBRAIC(:,48).*STATES(:,25)) - ( ALGEBRAIC(:,40).*STATES(:,20)+ ALGEBRAIC(:,46).*STATES(:,20));
ALGEBRAIC(:,50) = ALGEBRAIC(:,46)./1000.00;
ALGEBRAIC(:,52) = ALGEBRAIC(:,42);
RATES(:,25) = ( ALGEBRAIC(:,46).*STATES(:,20)+ ALGEBRAIC(:,44).*STATES(:,21)+ ALGEBRAIC(:,52).*STATES(:,23)+ ALGEBRAIC(:,27).*STATES(:,26)) - ( ALGEBRAIC(:,48).*STATES(:,25)+ ALGEBRAIC(:,42).*STATES(:,25)+ ALGEBRAIC(:,50).*STATES(:,25)+ ALGEBRAIC(:,38).*STATES(:,25));
i_pCa = ( i_pCa_max.*(Ca_i .^ 2.00000))./((Km_pCa.^ 2.00000)+(Ca_i .^ 2.00000));
i_NaCa =  (( (( (( k_NaCa.*1.00000)./((K_mNa .^ 3.00000)+(Na_o .^ 3.00000))).*1.00000)./(K_mCa+Ca_o)).*1.00000)./(1.00000+ k_sat.*(exp((( (eta - 1.00000).*V_m.*F)./( R.*T)))))).*( (exp((( eta.*V_m.*F)./( R.*T)))).*(Na_i .^ 3.00000).*Ca_o -  (exp((( (eta - 1.00000).*V_m.*F)./( R.*T)))).*(Na_o .^ 3.00000).*Ca_i);
E_CaN =  (( R.*T)./( 2.00000.*F)).*(log((Ca_o./Ca_i)));
i_CaB =  G_CaB.*(V_m - E_CaN);
B_i = (1.00000+( CMDN_tot.*Km_CMDN)./((Km_CMDN+Ca_i) .^ 2.00000)) .^ ( - 1.00000);
J_trpn = ( CONSTANTS(:,22).*Ca_i.*(CONSTANTS(:,33) - STATES(:,8))+ CONSTANTS(:,24).*Ca_i.*(CONSTANTS(:,32) - STATES(:,7))) - ( CONSTANTS(:,23).*STATES(:,8)+ CONSTANTS(:,25).*STATES(:,7));
Ca_i_prime =  B_i.*((J_leak+J_xfer) - (J_up+J_trpn+( ((i_CaB+i_pCa) -  2.00000.*i_NaCa).*A_cap.*C_m)./( 2.00000.*vol_myo.*F)));
ALGEBRAIC(:,54) = ALGEBRAIC(:,46)./95000.0;
ALGEBRAIC(:,56) = ALGEBRAIC(:,42)./50.0000;
RATES(:,23) = ( ALGEBRAIC(:,50).*STATES(:,25)+ ALGEBRAIC(:,56).*STATES(:,24)) - ( ALGEBRAIC(:,52).*STATES(:,23)+ ALGEBRAIC(:,54).*STATES(:,23));
RATES(:,24) =  ALGEBRAIC(:,54).*STATES(:,23) -  ALGEBRAIC(:,56).*STATES(:,24);
ALGEBRAIC(:,57) =  (( R.*T)./F).*(log((( 0.900000.*Na_o+ 0.100000.*CONSTANTS(:,7))./( 0.900000.*Na_i+ 0.100000.*STATES(:,28)))));
ALGEBRAIC(:,58) =  CONSTANTS(:,56).*STATES(:,20).*(V_m - ALGEBRAIC(:,57));
ALGEBRAIC(:,59) =  CONSTANTS(:,57).*(V_m - ALGEBRAIC(:,57));
ALGEBRAIC(:,68) = 1.00000./(1.00000+ 0.124500.*(exp(((  - 0.100000.*V_m.*F)./( R.*T))))+ 0.0365000.*CONSTANTS(:,72).*(exp(((  - V_m.*F)./( R.*T)))));
ALGEBRAIC(:,69) = ( (( CONSTANTS(:,66).*ALGEBRAIC(:,68).*1.00000)./(1.00000+((CONSTANTS(:,67)./Na_i) .^ 1.50000))).*CONSTANTS(:,7))./(CONSTANTS(:,7)+CONSTANTS(:,68));
RATES(:,19) = (  - (ALGEBRAIC(:,58)+ALGEBRAIC(:,59)+ 3.00000.*ALGEBRAIC(:,69)+ 3.00000.*i_NaCa).*A_cap.*C_m)./( vol_myo.*F);
ALGEBRAIC(:,60) =  (( R.*T)./F).*(log((CONSTANTS(:,7)./STATES(:,28))));
ALGEBRAIC(:,61) =  CONSTANTS(:,58).*(STATES(:,29) .^ 3.00000).*STATES(:,30).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,62) =  CONSTANTS(:,59).*STATES(:,31).*STATES(:,32).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,63) = ( (( 0.293800.*CONSTANTS(:,7))./(CONSTANTS(:,7)+210.000)).*(V_m - ALGEBRAIC(:,60)))./(1.00000+(exp(( 0.0896000.*(V_m - ALGEBRAIC(:,60))))));
ALGEBRAIC(:,64) =  CONSTANTS(:,60).*(STATES(:,33) .^ 2.00000).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,65) =  CONSTANTS(:,61).*STATES(:,34).*STATES(:,35).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,66) =  CONSTANTS(:,62).*STATES(:,36).*STATES(:,37).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,67) =  CONSTANTS(:,63).*STATES(:,38).*(V_m -  (( R.*T)./F).*(log((( 0.980000.*CONSTANTS(:,7)+ 0.0200000.*Na_o)./( 0.980000.*STATES(:,28)+ 0.0200000.*Na_i)))));
RATES(:,28) = (  - ((ALGEBRAIC(:,61)+ALGEBRAIC(:,62)+ALGEBRAIC(:,63)+ALGEBRAIC(:,64)+ALGEBRAIC(:,66)+ALGEBRAIC(:,65)+ALGEBRAIC(:,67)) -  2.00000.*ALGEBRAIC(:,69)).*A_cap.*C_m)./( vol_myo.*F);
i_stim = piecewise({t>=stim_start&t<=stim_end&(t - stim_start) -  (floor(((t - stim_start)./stim_period))).*stim_period<=stim_duration, stim_amp }, 0.00000);
ALGEBRAIC(:,70) = 0.200000./(1.00000+(exp(( - (V_m - 46.7000)./7.80000))));
ALGEBRAIC(:,71) =  (( CONSTANTS(:,69).*ALGEBRAIC(:,70).*Ca_i)./(Ca_i+CONSTANTS(:,71))).*(V_m - CONSTANTS(:,70));
V_m_prime =  - (i_CaL+i_pCa+i_NaCa+i_CaB+ALGEBRAIC(:,58)+ALGEBRAIC(:,59)+ALGEBRAIC(:,69)+ALGEBRAIC(:,61)+ALGEBRAIC(:,62)+ALGEBRAIC(:,63)+ALGEBRAIC(:,64)+ALGEBRAIC(:,65)+ALGEBRAIC(:,66)+ALGEBRAIC(:,67)+ALGEBRAIC(:,71)+i_stim);
RATES = RATES';
end

% Calculate algebraic variables
function ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, t)
P_C1 = 1.00000 - (P_C2+P_O1+P_O2);
ALGEBRAIC(:,5) =  0.180640.*(exp(( 0.0357700.*(V_m+30.0000))));
ALGEBRAIC(:,15) =  0.395600.*(exp((  - 0.0623700.*(V_m+30.0000))));
ALGEBRAIC(:,6) = ( 0.000152000.*(exp(( - (V_m+13.5000)./7.00000))))./( 0.00670830.*(exp(( - (V_m+33.5000)./7.00000)))+1.00000);
ALGEBRAIC(:,16) = ( 0.000950000.*(exp(((V_m+33.5000)./7.00000))))./( 0.0513350.*(exp(((V_m+33.5000)./7.00000)))+1.00000);
ALGEBRAIC(:,7) = 1.00000./(1.00000+(exp(( - (V_m+22.5000)./7.70000))));
ALGEBRAIC(:,17) =  0.493000.*(exp((  - 0.0629000.*V_m)))+2.05800;
ALGEBRAIC(:,8) = 1.00000./(1.00000+(exp(((V_m+45.2000)./5.70000))));
ALGEBRAIC(:,18) = 270.000+1050.00./(1.00000+(exp(((V_m+45.2000)./5.70000))));
ALGEBRAIC(:,9) = ( 4.81333e-06.*(V_m+26.5000))./(1.00000 - (exp((  - 0.128000.*(V_m+26.5000)))));
ALGEBRAIC(:,19) =  9.53333e-05.*(exp((  - 0.0380000.*(V_m+26.5000))));
ALGEBRAIC(:,20) =  0.493000.*(exp((  - 0.0629000.*V_m)))+2.05800;
ALGEBRAIC(:,21) = 1200.00 - 170.000./(1.00000+(exp(((V_m+45.2000)./5.70000))));
ALGEBRAIC(:,22) =  39.3000.*(exp((  - 0.0862000.*V_m)))+13.1700;
ALGEBRAIC(:,11) =  0.0137330.*(exp(( 0.0381980.*V_m)));
ALGEBRAIC(:,24) =  6.89000e-05.*(exp((  - 0.0417800.*V_m)));
C1 = 1.00000 - (O+C2+C3+C4+I1+I2+I3);
alpha = ( 0.400000.*(exp(((V_m+12.0000)./10.0000))).*((1.00000+ 0.700000.*(exp(( - ((V_m+40.0000) .^ 2.00000)./10.0000)))) -  0.750000.*(exp(( - ((V_m+20.0000) .^ 2.00000)./400.000)))))./(1.00000+ 0.120000.*(exp(((V_m+12.0000)./10.0000))));
beta =  0.0500000.*(exp(( - (V_m+12.0000)./13.0000)));
ALGEBRAIC(:,10) = 1.00000 - (STATES(:,39)+STATES(:,40)+STATES(:,38)+STATES(:,41));
ALGEBRAIC(:,23) =  0.0223480.*(exp(( 0.0117600.*V_m)));
ALGEBRAIC(:,28) =  0.0470020.*(exp((  - 0.0631000.*V_m)));
ALGEBRAIC(:,29) =  0.0908210.*(exp(( 0.0233910.*(V_m+5.00000))));
ALGEBRAIC(:,33) =  0.00649700.*(exp((  - 0.0326800.*(V_m+5.00000))));
gamma = ( k_pc_max.*Ca_ss)./(k_pc_half+Ca_ss);
K_pcf =  13.0000.*(1.00000 - (exp(( - ((V_m+14.5000) .^ 2.00000)./100.000))));
B_JSR = (1.00000+( CSQN_tot.*Km_CSQN)./((Km_CSQN+Ca_JSR) .^ 2.00000)) .^ ( - 1.00000);
J_rel =  v1.*(P_O1+P_O2).*(Ca_JSR - Ca_ss).*P_RyR;
J_tr = (Ca_NSR - Ca_JSR)./tau_tr;
J_leak =  v2.*(Ca_NSR - Ca_i);
J_up = ( v3.*(Ca_i .^ 2.00000))./((K_m_up .^ 2.00000)+(Ca_i .^ 2.00000));
ALGEBRAIC(:,4) = 1.00000 - (STATES(:,20)+STATES(:,21)+STATES(:,22)+STATES(:,25)+STATES(:,23)+STATES(:,24)+STATES(:,26)+STATES(:,27));
ALGEBRAIC(:,14) = 3.80200./( 0.102700.*(exp(( - (V_m+2.50000)./17.0000)))+ 0.200000.*(exp(( - (V_m+2.50000)./150.000))));
ALGEBRAIC(:,36) =  0.191700.*(exp(( - (V_m+2.50000)./20.3000)));
ALGEBRAIC(:,27) = 3.80200./( 0.102700.*(exp(( - (V_m+2.50000)./15.0000)))+ 0.230000.*(exp(( - (V_m+2.50000)./150.000))));
ALGEBRAIC(:,38) =  0.200000.*(exp(( - (V_m - 2.50000)./20.3000)));
ALGEBRAIC(:,42) =  7.00000e-07.*(exp(( - (V_m+7.00000)./7.70000)));
ALGEBRAIC(:,44) = 0.00840000+ 2.00000e-05.*(V_m+7.00000);
ALGEBRAIC(:,32) = 3.80200./( 0.102700.*(exp(( - (V_m+2.50000)./12.0000)))+ 0.250000.*(exp(( - (V_m+2.50000)./150.000))));
ALGEBRAIC(:,40) =  0.220000.*(exp(( - (V_m - 7.50000)./20.3000)));
i_CaL =  G_CaL.*O.*(V_m - E_CaL);
B_ss = (1.00000+( CMDN_tot.*Km_CMDN)./((Km_CMDN+Ca_ss) .^ 2.00000)) .^ ( - 1.00000);
J_xfer = (Ca_ss - Ca_i)./tau_xfer;
ALGEBRAIC(:,46) = 1.00000./( 0.188495.*(exp(( - (V_m+7.00000)./16.6000)))+0.393956);
ALGEBRAIC(:,48) = ( ALGEBRAIC(:,32).*ALGEBRAIC(:,46).*ALGEBRAIC(:,42))./( ALGEBRAIC(:,40).*ALGEBRAIC(:,44));
ALGEBRAIC(:,50) = ALGEBRAIC(:,46)./1000.00;
ALGEBRAIC(:,52) = ALGEBRAIC(:,42);
i_pCa = ( i_pCa_max.*(Ca_i .^ 2.00000))./((Km_pCa.^ 2.00000)+(Ca_i .^ 2.00000));
i_NaCa =  (( (( (( k_NaCa.*1.00000)./((K_mNa .^ 3.00000)+(Na_o .^ 3.00000))).*1.00000)./(K_mCa+Ca_o)).*1.00000)./(1.00000+ k_sat.*(exp((( (eta - 1.00000).*V_m.*F)./( R.*T)))))).*( (exp((( eta.*V_m.*F)./( R.*T)))).*(Na_i .^ 3.00000).*Ca_o -  (exp((( (eta - 1.00000).*V_m.*F)./( R.*T)))).*(Na_o .^ 3.00000).*Ca_i);
E_CaN =  (( R.*T)./( 2.00000.*F)).*(log((Ca_o./Ca_i)));
i_CaB =  G_CaB.*(V_m - E_CaN);
B_i = (1.00000+( CMDN_tot.*Km_CMDN)./((Km_CMDN+Ca_i) .^ 2.00000)) .^ ( - 1.00000);
J_trpn = ( CONSTANTS(:,22).*Ca_i.*(CONSTANTS(:,33) - STATES(:,8))+ CONSTANTS(:,24).*Ca_i.*(CONSTANTS(:,32) - STATES(:,7))) - ( CONSTANTS(:,23).*STATES(:,8)+ CONSTANTS(:,25).*STATES(:,7));
ALGEBRAIC(:,54) = ALGEBRAIC(:,46)./95000.0;
ALGEBRAIC(:,56) = ALGEBRAIC(:,42)./50.0000;
ALGEBRAIC(:,57) =  (( R.*T)./F).*(log((( 0.900000.*Na_o+ 0.100000.*CONSTANTS(:,7))./( 0.900000.*Na_i+ 0.100000.*STATES(:,28)))));
ALGEBRAIC(:,58) =  CONSTANTS(:,56).*STATES(:,20).*(V_m - ALGEBRAIC(:,57));
ALGEBRAIC(:,59) =  CONSTANTS(:,57).*(V_m - ALGEBRAIC(:,57));
ALGEBRAIC(:,68) = 1.00000./(1.00000+ 0.124500.*(exp(((  - 0.100000.*V_m.*F)./( R.*T))))+ 0.0365000.*CONSTANTS(:,72).*(exp(((  - V_m.*F)./( R.*T)))));
ALGEBRAIC(:,69) = ( (( CONSTANTS(:,66).*ALGEBRAIC(:,68).*1.00000)./(1.00000+((CONSTANTS(:,67)./Na_i) .^ 1.50000))).*CONSTANTS(:,7))./(CONSTANTS(:,7)+CONSTANTS(:,68));
ALGEBRAIC(:,60) =  (( R.*T)./F).*(log((CONSTANTS(:,7)./STATES(:,28))));
ALGEBRAIC(:,61) =  CONSTANTS(:,58).*(STATES(:,29) .^ 3.00000).*STATES(:,30).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,62) =  CONSTANTS(:,59).*STATES(:,31).*STATES(:,32).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,63) = ( (( 0.293800.*CONSTANTS(:,7))./(CONSTANTS(:,7)+210.000)).*(V_m - ALGEBRAIC(:,60)))./(1.00000+(exp(( 0.0896000.*(V_m - ALGEBRAIC(:,60))))));
ALGEBRAIC(:,64) =  CONSTANTS(:,60).*(STATES(:,33) .^ 2.00000).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,65) =  CONSTANTS(:,61).*STATES(:,34).*STATES(:,35).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,66) =  CONSTANTS(:,62).*STATES(:,36).*STATES(:,37).*(V_m - ALGEBRAIC(:,60));
ALGEBRAIC(:,67) =  CONSTANTS(:,63).*STATES(:,38).*(V_m -  (( R.*T)./F).*(log((( 0.980000.*CONSTANTS(:,7)+ 0.0200000.*Na_o)./( 0.980000.*STATES(:,28)+ 0.0200000.*Na_i)))));
i_stim = piecewise({t>=stim_start&t<=stim_end&(t - stim_start) -  (floor(((t - stim_start)./stim_period))).*stim_period<=stim_duration, stim_amp }, 0.00000);
ALGEBRAIC(:,70) = 0.200000./(1.00000+(exp(( - (V_m - 46.7000)./7.80000))));
ALGEBRAIC(:,71) =  (( CONSTANTS(:,69).*ALGEBRAIC(:,70).*Ca_i)./(Ca_i+CONSTANTS(:,71))).*(V_m - CONSTANTS(:,70));
end

% Compute result of a piecewise function
function x = piecewise(cases, default)
set = [0];
for i = 1:2:length(cases)
    if (length(cases{i+1}) == 1)
        x(cases{i} & ~set,:) = cases{i+1};
    else
        x(cases{i} & ~set,:) = cases{i+1}(cases{i} & ~set);
    end
    set = set | cases{i};
    if(set), break, end
end
if (length(default) == 1)
    x(~set,:) = default;
else
    x(~set,:) = default(~set);
end
end

% Pad out or shorten strings to a set length
function strout = strpad(strin)
req_length = 160;
insize = size(strin,2);
if insize > req_length
    strout = strin(1:req_length);
else
    strout = [strin, blanks(req_length - insize)];
end
end


